function [ Z ] = getKernelGradient( p, q )
%EVALUATEQUADRATUREPOINT Summary of this function goes here
%   Detailed explanation goes here

% shiftAlpha = 1e-3;

global shiftAlpha shearModulus poissonRatio;

r = p - q;

rn = max(norm(r), shiftAlpha);

rg = r / rn;

M = [rg(1)*rg(1), rg(1)*rg(2), rg(1)*rg(3); rg(2)*rg(1), rg(2)*rg(2), rg(2)*rg(3); rg(3)*rg(1), rg(3)*rg(2), rg(3)*rg(3)];
I = diag(ones(3, 1));

DrM = getDrM(r, rn, M);
DrI = getDrI(r, rn, I);

Z = ((3 - 4 * poissonRatio) * DrI + DrM) / (16 * pi * shearModulus * (1 - poissonRatio));

% Z = DrM;

% Z = (DrI + DrM) / (16 * pi);

end

